%% DERIVATIVES 

Param = readmatrix('Sigmastar_1sector.csv');
pm.sigma = Param(3);
pm.se = Param(4);
pm.theta = Param(5);
pm.phi = Param(6);
pm.bc = Param(7);
pm.bp = Param(8);
pm.E_i = Param(9);
pm.E_r = Param(10);
pm.price = Param(11);
pm.wage = Param(12);

fprintf('Derivatives \n');
adjfactor = Inf;
base = pm;

%Baseline
fprintf('Baseline \n');
i = 1;
[p,w,~,ZstarMat, LMat,~,NMat,~,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
zgrid = linspace(1,2,1000)'.*ZstarMat(:,2);
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

writematrix([data.ind data.NIC ...
             mean(cr)./(1+mean(cr)) NMat LMat./NMat VMat' coeff(1)] , ...
             ['Derivatives_1s_' num2str(i) '.csv'],'delimiter',',');


%Sigma x
fprintf('Sigma_x \n');
i = i+1;
pm.sigma = 1.01.*base.sigma;

[p,w,~,ZstarMat, LMat,~,NMat,~,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
zgrid = linspace(1,2,1000)'.*ZstarMat(:,2);
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

writematrix([data.ind data.NIC ...
             mean(cr)./(1+mean(cr)) NMat LMat./NMat VMat' coeff(1)] , ...
             ['Derivatives_1s_' num2str(i) '.csv'],'delimiter',',');

%Sigma e
fprintf('Sigma_e \n');
i = i+1;
pm = base;
pm.se = 1.01.*base.se;

[p,w,~,ZstarMat, LMat,~,NMat,~,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
zgrid = linspace(1,2,1000)'.*ZstarMat(:,2);
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

writematrix([data.ind data.NIC ...
             mean(cr)./(1+mean(cr)) NMat LMat./NMat VMat' coeff(1)] , ...
             ['Derivatives_1s_' num2str(i) '.csv'],'delimiter',',');

%Theta 
fprintf('Theta \n');
i = i+1;
pm = base;
pm.theta = 1.01.*base.theta;

[p,w,~,ZstarMat, LMat,~,NMat,~,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
zgrid = linspace(1,2,1000)'.*ZstarMat(:,2);
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

writematrix([data.ind data.NIC ...
             mean(cr)./(1+mean(cr)) NMat LMat./NMat VMat' coeff(1)] , ...
             ['Derivatives_1s_' num2str(i) '.csv'],'delimiter',',');

%Phi 
fprintf('Phi \n');
i = i+1;
pm = base;
pm.phi= 1.01.*base.phi;

[p,w,~,ZstarMat, LMat,~,NMat,~,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
zgrid = linspace(1,2,1000)'.*ZstarMat(:,2);
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

writematrix([data.ind data.NIC ...
             mean(cr)./(1+mean(cr)) NMat LMat./NMat VMat' coeff(1)] , ...
             ['Derivatives_1s_' num2str(i) '.csv'],'delimiter',',');


%bc
fprintf('Bc \n');
i = i+1;
pm = base;
pm.bc = 1.01.*base.bc;

[p,w,~,ZstarMat, LMat,~,NMat,~,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
zgrid = linspace(1,2,1000)'.*ZstarMat(:,2);
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

writematrix([data.ind data.NIC ...
             mean(cr)./(1+mean(cr)) NMat LMat./NMat VMat' coeff(1)] , ...
             ['Derivatives_1s_' num2str(i) '.csv'],'delimiter',',');

%bp
fprintf('Bp \n');
i = i+1;
pm = base;
pm.bp = 1.01.*base.bp;

[p,w,~,ZstarMat, LMat,~,NMat,~,VMat] = Counterfactuals(pm,pm.bc,pm.bp);
zgrid = linspace(1,2,1000)'.*ZstarMat(:,2);
lnL = log(l_f(pm,zgrid,1,p,w,pm.t,1));
cr = (pm.bc./(zgrid.^(pm.phi).*pm.bp)).^(-pm.nu);
cr = 1 - (1./(1+cr));
coeff = polyfit(lnL,cr,2);

writematrix([data.ind data.NIC ...
             mean(cr)./(1+mean(cr)) NMat LMat./NMat VMat' coeff(1)] , ...
             ['Derivatives_1s_' num2str(i) '.csv'],'delimiter',',');

fprintf('\n \n Done: Derivatives \n \n');
